
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include<sstream>

#include"library.h"

#define OmegaM 0.286
#define OmegaA 0.714
#define c 3E10
const double H0 = 69.6E-19 / 3.08568 ;
const double  Lp =  2.74E42 * 1E23;
#define vc 6E8

const double integ_eps = 1E-3;

double g(double zz){
	double y = pow(1+zz, 3) * OmegaM + OmegaA;
	return 1.0/sqrt(y);
}

double DL( double z ){
	double integ = SelfAdapt_Simpson( 0, z, g, integ_eps );
	return c / H0 * (1+z) * integ;
}

double S( double z ){
	double y = DL( z );
	return Lp / 4 / M_PI / y / y * pow(1+z, 1.0/3) / vc;
}

int main(){
	
	string mode;
	cout<<" Mode Options: \n";
	cout<<" (DL) calculate DL(z), with given z \n";
	cout<<" (S) calculate S(z), with given z \n";
	cout<<" (dSdz) calculate dSdz, wigh given z \n";
	cout<<" (plot) plot log10(dSdz) ~ log10(1+z) \n";
	string line;
	stringstream strm; strm.clear(); getline(cin, line); strm.str(line);
	strm >> mode;
	cout<<"mode = "<<mode<<endl;

	double integ_eps = 1E-3, z;
	double h = 0.01;

	if( mode == "DL" ){
		cout<<" Please input value of z = ";
		strm.clear(); getline(cin, line); strm.str(line); strm >> z;
		cout<<" integ precision = "<<integ_eps<<", z = "<< z <<", DL(z) = "<<DL(z)<<endl;
	}
	else if( mode == "S" ){
		cout<<" Please input value of z = ";
		strm.clear(); getline(cin, line); strm.str(line); strm >> z;
		cout<<" integ precision = "<<integ_eps<<", z = "<< z <<", S(z) = "<< S(z)<<endl;
	}
	else if( mode == "dSdz" ){
		cout<<" Please input value of z = ";
		strm.clear(); getline(cin, line); strm.str(line); strm >> z;
		cout<<" integ precision = "<<integ_eps<<", z = "<< z <<", dSdz(z) = "<< Richardson_deriv(z, h, 3, S)<<endl;
	}
	else if( mode == "plot" ){
		h = 0.001;
		ofstream fout;
		fout.open("dSdz.txt");
		for(double z = 2*h; z<1; z += h )
			fout<<log10(1+z)<<"  "<<log10( - Richardson_deriv( z, h, 3, S ) )<<endl;
		for(double z = 1; z<10; z += 0.001 )
			fout<<log10(1+z)<<"  "<<log10( - Richardson_deriv( z, h, 3, S ) )<<endl;
		fout.close();
		cout<<" log10( - dSdz ) ~ log10( 1+z ) data is written into dSdz.txt \n";
	}

	return 0;
}
